Classical potential describes martensitic phase transformations between the a, (3 and 

oj titanium phases. 
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A description of the martensitic transformations between the a, (3 and uj phases of titanium 
that includes nucleation and growth requires an accurate classical potential. Optimization of the 
parameters of a modified embedded atom potential to a database of density-functional calculations 
yields an accurate and transferable potential as verified by comparison to experimental and density 
functional data for phonons, surface and stacking fault energies and energy barriers for homogeneous 
martensitic transformations. Molecular dynamics simulations map out the pressure-temperature 
phase diagram of titanium. For this potential the martensitic phase transformation between a and 
appears at ambient pressure and 1200 K, between a and uj at ambient conditions, between f3 and 
uj at 1200 K and pressures above 8 GPa, and the triple point occurs at 8GPa and 1200 K. Molecular 
dynamics explorations of the dynamics of the martensitic a — uj transformation show a fast-moving 
interface with a low interfacial energy of 30 meV/A 2 . The potential is applicable to the study of 
defects and phase transformations of Ti. 
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I. INTRODUCTION 

Martensitic phase transitions control systems rang- 
ing from shape memory alloys^ to steels^ to planetary 
cores.'' They are diffusionless structural transformations 
proceeding near the speed of sound^ Martensitic trans- 
formations frequently appear in alloy design as a way to 
improve materials properties, but their occurrence can 
also limit materials performance. 

Titanium's great technological importance^ makes it 
an ideal example for the development of physics-based 
predictive methods for materials problems. Titanium 
displays several phases as a function of pressure and tem- 
perature. At ambient conditions titanium stabilizes in 
the hexagonal close-packed (hep) a phase. At ambient 
pressure and temperatures above 1155 K the a phase 
transforms to the high-temperature body-centered cubic 
(bec) (3 phase. Under pressure the a phase transforms 
into the hexagonal u> phased The high-pressure uj phase 
consist of a three atom hexagonal structure equivalent 
to AIB2 with Ti occupying sites of alternating triangu- 
lar and honeycomb layers. The crystal structure of Ti at 
K has not been determined experimentally. Extrapola- 
tion of the a-uj phase boundary^ in Ti indicates uj as the 
ground state phase. Free energy calculations of a and 
uj within the quasiharmonic approximation using a TB 
model show phonon entropy stabilizing the a phase at 
ambient temperature^ Ti displays martensitic transfor- 
mations between its a (hep), (3 (bec) and uj (hexagonal) 
phases. Two-phase a/ (3 alloys make up many industrial 
titanium alloys such as Ti-6A1-4V because the presence 
of (3 phase in the a matrix improves strength and creep 



resistance. Titanium transforms from a to brittle uj un- 
der pressure creating serious technological problems for 
(3 stabilized titanium alloys. Impurities greatly affect the 
a to uj transformation; for example, as little as 1 at.% 
oxygen in commercial Ti alloys suppresses it^^^ 

To understand these transformations in Ti and its al- 
loys we begin with the study of the phase transformations 
in pure titanium. Our approach to a theoretical under- 
standing of these transformation involves three steps: (1) 
find the homogeneous atomistic pathway of the marten- 
sitic transformation in pure titanium, (2) use this path- 
way to estimate the effect of impurities, and (3) deter- 
mine the heterogeneous nucleation and dynamics of the 
martensitic phase transformation. 

The homogeneous transformation pathway for all three 
martensitic transformations in Ti is known. Burgers de- 
scribed the a to (3 transformation in Zr j 1 ^ the same mech- 
anism occurs in Ti. The (3 to uj transformation occurs 
via plane collapse along the [111] direction corresponding 
to the longitudinal §[111] phonom 12 : 13 ' 14 More recently 
Trinkle et al. determined the homogeneous pathway of 
the a to uj transformation ! 15 ' 16 A systematic approach 
generated all possible pathways that were then succes- 
sively pruned by energy estimates using elastic theory, 
tight-binding (TB) methods and density-functional the- 
ory (DFT). 

The speed of the diffusionless martensitic transforma- 
tion traps dilute impurities, providing candidate path- 
ways for alloyed materials. Hennig et al. determined the 
effect of interstitial and substitutional impurities on the 
a to w transformation. 10 DFT nudgcd-elastic band re- 
finements yield the change in both the relative stability 
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of and the energy barrier between the phases due to im- 
purities. The resulting microscopic picture explains the 
suppression of the a to ui transformation in commercial 
Ti alloys. 

The final step involves studying the full atomistic dy- 
namics of the nucleation and growth of the martensitic 
phase; this requires molecular dynamics simulations of 
large systems. For the required system sizes an accurate 
quantum mechanical treatment by DFT or TB methods 
becomes too computationally demanding. Such simula- 
tions call for a classical potential to allow the treatment 
of appropriate length and time scales for nucleation and 
growth of martensites. 

In this paper we develop a classical potential of the 
modified embedded atom method^ (MEAM) type for Ti. 
The potential accurately describes the stability of the a, 
j3 and uj phases and is applied to study the dynamics of 
the martensitic atow transformation and the interfacial 
energy between the a and ui phase. Section HT1 describes 
the calculations for the DFT database, the functional 
form of the MEAM potential and the optimization of the 
potential parameters to the DFT database. The accu- 
racy of the potential is tested by comparing phonon spec- 
tra, surface and stacking fault energies as well as energy 
barriers for homogeneous martensitic transformations to 
DFT, TB and experimental results. In Section IIIII we 
apply the potential to study the phase diagram of Ti 
and the martensitic phase transformations between the 
phases. We estimate the interfacial energy between a and 
u> and show that the classical MEAM potential accurately 
describes the stability range of the three Ti phases and 
the phase transformations between them. 



II. OPTIMIZATION OF THE CLASSICAL 
POTENTIAL TO DENSITY FUNCTIONAL 
DATABASE 



To describe the interactions between the Ti atoms and 
to enable large-scale molecular dynamics simulations we 
develop a classical potential. The modified embedded 
atom method provides the form of the potential™ with 
potential parameters optimized to a database of DFT 
calculations. A second DFT database provides testing 
data for the potential. The optimization of the model 
proceeds iteratively. Systematically adding DFT results 
to the fitting and testing databases improves the accu- 
racy and extends the applicability of the model. This 
enables the development of a potential that accurately 
reproduces the properties of all three Ti phases relevant 
for the description of the martensitic phase transitions. 
Available experimental data confirms the accuracy of the 
resulting Ti MEAM potential. 



A. Density functional calculations 



The DFT calculations are performed with VASP ) 18 i 19 a 
density functional code using a plane-wave basis and ul- 
trasoft Vanderbilt type pseudopotent ials j 20 i 2 1 The gen- 
eralized gradient approximation (GGA) of Perdew and 
Wang is used^ 2 - A plane-wave kinetic energy cutoff of 
400 eV ensures energy convergence to 0.3 meV/atom. 
The fc-point meshes for the different structures are cho- 
sen to guarantee an energy accuracy of 1 meV/atom. We 
treat the Ti 3p states as valence states in addition to the 
usual 4s and 3d states to provide an accurate treatment 
of the interaction at close interatomic distances. 

The DFT database for the fitting of the potential pa- 
rameters consists of energies, defects, forces and elastic 
constants for a variety of Ti phases as well as energies of 
configurations along the TAO-1 transformation pathway 
from a to uji 15 ^ 16 Relaxations determine the ground state 
energies and lattice parameters of the a, j3, uj, fee, A15 
and simple hexagonal phases. The volume dependence 
of the energy for a, f3, uj, fee and A15 and the elastic 
constants of the a, (3 and u> phases at their equilibrium 
volumes are calculated. Snapshots of short DFT molec- 
ular dynamics simulations for a, [3 and w at 800 K at 
their respective ground state volumes with and without 
a vacancy defect provide force-matching data: 23 Defect 
structures and formation energies are determined by re- 
laxations for single vacancies and single interstitial atoms 
in a 96-atom (4x4x3) supercell for a and a 108-atom 
(3x3x4) supercell for u with a 2 x 2 x 2 fc-point sampling 
grid. This results in a 1 at.% defect concentration. The 
atom positions are relaxed until the atomic-level forces 
are smaller than 20 meV/A. 

The DFT database for the test of the accuracy of the 
potential consists of additional interstitial defects in a 
and uj, phonon spectra for a, (3 and u>, surface energies 
for a and uj, and the I2 stacking fault in a. The phonon 
calculation employ the direct force method and super- 
cells of 150 atoms (5x5x3) for a, 125 atoms (5x5x5) 
for (3, and 135 atoms (3x3x5) for uj. The surfaces are 
constructed by separating the crystal along a high-index 
plane. The surface energies result from relaxing a peri- 
odic stacking of an 18 A to 20 A thick slab of rectangular 
hep cells with a 10 A vacuum region and a perfect bulk 
cell with the same cell vectors. A /c-point mesh equiv- 
alent to 13 x 13 x 1 is used for both the bulk and the 
slab calculation. Relaxations of a 1 x 1 x 5 supercell of 
a with and without a single I2 stacking fault determine 
the stacking fault energy. 

Comparison to available experimental data for phonon 
spectra, surface energies, the stacking fault in a, and the 
p-T phase diagram of Ti further confirm the accuracy 
of the potential as do approximate energy barriers for 
homogeneous martensitic transformations in TB. 
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FIG. 1: The five cubic splines of the MEAM potential. The lines denote the spline interpolation between the spline knots 
represented by the solid circles. The splines are linearly extrapolated beyond the last spline knot. 



B. Modified embedded atom potential 

The MEAM formalism was originally developed by 
Basket as an extension of the embedded-atom method. 
The original MEAM includes an angular-dependent elec- 
tron density to model the effects of bond bending; a series 
of four terms with s, p, d and / character describe the an- 
gular densities. The original MEAM potential has been 
applied to a variety of systems ranging from the semi- 
conductors Si2ii2£ and Ge2& to bcc and fee metals2£>2£ 
to several binary alloy system s 26 ' 29 and recently to hep 
Ti and Zr£2- While the Ti and Zr potentials accurately 
reproduce properties of the hep phase, they do not de- 
scribe the bcc and w phases^ Here we aim to develop a 
potential that accurately describes all three phases and 
the transformations between them. 

More recently, Lenosky et al. modified the original 
MEAM potential by using cubic splines for the functional 
fornijii This removes the constraint of fixed angular char- 
acter and allows for additional flexibility of the potential. 
Furthermore, the use of splines reduces the cost of evalua- 
tion over the original functional form providing increased 
computational efficiency In practice, the evaluation of 
the spline-based MEAM potentials is only about a factor 
of two slower than that of EAM potentials. The spline- 
based MEAM was successfully applied to study Si.— £^ 
This success of the spline-based MEAM, its improved 
flexibility and its higher computational efficiency moti- 
vate our use of this functional form here. The MEAM 
potential is implemented into two freely-available large- 
scale parallel molecular dynamics codesi 33 ' 34 

The MEAM potential used in this work separates the 
energy into two parts^i 

ij i 

with the density at atom i 

Pi =J2p( r *j) + E /( r u)/( r »fe)5(cos(% fe )), (2) 

3 jk 

where Ojik is the angle between atoms j, i and k cen- 
tered on atom i. The functional form contains as special 
cases the Stillinger- Weber— (U(x) = x and p = 0) and 
the embedded-atom (EAM) (/ = or g = 0) potentials. 



TABLE I: Parameters specifying the five cubic splines that 
comprise the MEAM potential. The first part of the table 
lists the number of knots N for each spline and the range of 
the spline variables i m i n and f max . The middle part of the 
table gives the values at equally spaced spline knots defined 
by ti — t m i n + i(t max — tmi n )/N where N is the number of 
spline knots. Finally, the derivatives of the splines at their 
endpoints are listed in the last part of the table. 
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r[A] 


1.7427 


5.5000 


13 
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r[A] 


2.0558 


4.4100 


11 
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r[A] 


2.0558 


4.4100 


10 
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Ptot 


-55.1423 


-23.9383 


4 
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cos(0) 


-1.0000 


0.9284 


8 
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4>{n) [eV] 


p(ri) 


fin) 


U{ Pi ) [eV] 







3.7443 


1.7475 


-0.1485 


-0.29746 


0.0765 
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0.9108 


-5.8678 


1.6845 


-0.15449 


0.1416 


2 


0.3881 


-8.3376 


2.0113 


0.05099 


0.7579 
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-0.0188 


-5.8399 


1.1444 


0.57343 


0.6301 


4 


-0.2481 


-3.1141 


0.2862 




0.0905 


5 


-0.2645 


-1.7257 


-0.3459 




-0.3574 


6 


-0.2272 


-0.4429 


-0.6258 




-0.6529 


7 


-0.1293 


-0.1467 


-0.6120 




-6.0091 


8 


-0.0597 


-0.2096 


-0.3112 
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-0.0311 


-0.1442 


0.0000 






10 


-0.0139 


0.0000 
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-0.0032 
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0.0000 
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p'(r) 


fir) 


U'(p) 


g'ip) 




[eV/A] 


[A" 1 ] 


[A- 1 ] 


[eV] 







-20.0 


-1.0 


2.7733 


0.0078 


8.3364 


N 


0.0 


0.0 


0.0000 


0.1052 


-60.4025 



The five functions <5(r), U(p), p(r), f(r) and g(cos(8j) 
are represented by cubic splines ^ This allows for the 
necessary flexibility to accurately describe a complex sys- 
tem such as Ti and provides the computational efficiency 
required for large scale molecular dynamics simulations. 



C. MEAM Potential fit 

The spline parameters are optimized by a novel algo- 
rithm that involves an extensive parameter search. A 
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Volume [A 3 ] Pressure t GPa ] 

FIG. 2: (color online) Fitted energies of a, /3 and uj Ti as 
a function of volume, (a) The MEAM potential is fit to the 
energy of a, j3 and uj at several volumes. The DFT results 
appear as symbols and the fitted MEAM curves as lines, (b) 
The enthalpy difference between a and uj as a function of pres- 
sure for MEAM. The MEAM potential predicts a transition 
pressure p a -ui of about —5 GPa at zero temperature similar 
to DFT which predicts the uj phase to be the ground state.— 



TABLE II: Fitted energies relative to the a phase and lat- 
tice parameters for different Ti phases within the MEAM 
potential and DFT calculations (in parentheses). The co- 
hesive energy of the a. phase is 5.171 eV/atom in GGA and 
4.831 eV/atom in MEAM. 



Phase _E c [meV/atom] 


V [A 3 ] 


c/a 


a (0) 
uj -5 (-5) 
P 111 (108) 


17.40 (17.55) 
17.24 (17.28) 
17.51 (17.34) 


1.596 (1.583) 
0.611 (0.619) 


fee 39 (58) 
A15 54 (192) 
hexagonal 396 (353) 


17.83 (17.53) 
17.99 (17.49) 
17.25 (17.78) 


0.982 (0.999) 



detailed description of the algorithm will be published 
separately.— Figure [1] and Table UJ show the splines and 
spline parameters of the best MEAM potential. 

Predictions of the resulting potential confirm its ac- 
curacy and transferability. Table [Til compares the DFT 
energies and lattice parameters with the MEAM values 
for the experimentally observed a, and uj phases and 
the fee, A15 and simple hexagonal structures. The fit- 
ted MEAM values of the cohesive energy of the a phase 
of 4.831 eV and the lattice parameter of 2.931 A agree 
closely with the experimental values (4.844 eV, 2.951 A) 
and the DFT results (5.171 eV, 2.948 A). For the two 
low-energy structures a and uj the final MEAM potential 
is also fitted to the energy as a function of volume. For 
and fee the fit includes the equilibrium lattice constants 
and energies relative to hep. For the simple hexagonal 
structure the potential is fitted to the energy relative to 
hep for the structure with both lattice parameters fixed 
to the DFT values. 

Figure [^a) shows the result of the energy fit at differ- 
ent volumes for a, (3 and uj. The DFT and MEAM ener- 
gies agree to about 5 meV/atom with closer agreement 
near the energy minimum and slightly worse agreement 
at high compression. The MEAM potential reproduces 



TABLE III: Fitted elastic constants in units of GPa for a, P 
and uj Ti for the MEAM potential compared to the DFT fit- 
ting data and experiment. The MEAM potential accurately 
reproduces the DFT elastic constants. The experimental val- 
ues for a Ti are measured at 4 K (~R.ef.l38h and the P Ti elastic 
constants at 1238 K fR.ef.l39T). The deviation between the cal- 
culated and measured elastic constants for /3 Ti stems from 
the high temperature needed to stabilize the structure. For 
a and uj Ti, internal relaxations are neglected in the MEAM 
and DFT calculations. 
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MEAM 
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95 


58 
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GGA 
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82 


45 
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75 


Exp. 
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87 

— uj -phase 


51 
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MEAM 
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64 


GGA 
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81 

— P-phase 


54 


245 


54 


MEAM 


95 


111 


53 






GGA 


95 
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42 






Exp. 


134 


110 


36 







the uj phase as the ground state and places the a phase 
slightly higher in energy by 5 meV/atom in agreement 
with DFT calculations. 7 Figure [^b) shows the resulting 
enthalpy difference between a and uj. At zero tempera- 
ture the a phase transforms to uj at a pressure of —5 GPa 
in the MEAM potential. 

Table IIIII compares the elastic constants in MEAM 
with DFT results and experiment for the a, (3 and uj 
phases. For both methods the elastic constants are cal- 
culated neglecting internal relaxations in the a and uj 
phases. Accurate elastic constants are important for the 
correct description of the long-range strain fields around 
dislocations and other defect structures as well as in 
martensitic transformations. The low RMS deviation 
between MEAM and DFT elastic constants of 13% and 
maximum deviation of 29% demonstrates the quality of 
the fit and indicates the accuracy of the potential for the 
effects of strain in Ti. 

Table HVl compares the formation energies of point de- 
fects in the a and uj phases for the MEAM potential with 
DFT and several TB potentials. The defect relaxations 
are performed at fixed equilibrium volume for single de- 
fects in a 4 x 4 x 3 and 3x3x4 supercell of a and w, 
respectively. This corresponds to a defect concentration 
of 1%. The formation energies of the various interstitial 
atoms and vacancies in both phases agree well with the 
DFT results. In fact, the RMS errors of the energies are 
similar for the classical MEAM and the TB potentials. 
In addition the MEAM potential stabilizes the correct 
defects found in DFT calculations. The fitting data did 
not include the a dumbbell- [0001], the a divacancy-AB, 
or the uj hexahedral interstitial defects. The close agree- 
ment for these defects indicates the model's accuracy. 
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Reduced wave vector coordinate £, 

FIG. 3: (color online) The phonon spectra of a, f3 and u Ti for the classical MEAM potential (black solid lines) compared to 
DFT results (red dashed lines) and to experimental data for the a and j3 phases (blue circles) . 39 ' 42 For the a phase the MEAM 
potential accurately reproduces the low energy acoustic branches reflecting the elastic constants. For optical and high-energy 
acoustic branches the potential reproduces the experimental phonon frequencies within 25%, generally underestimating the 
experimental values. Similarly for the ui phase the acoustic branches of the MEAM potential closely match the DFT results 
with larger discrepancies of up to 20% for the optical modes. The j3 phase is mechanically unstable 3 -^ at low temperatures and 
transforms to a or uj. This instability is reflected in unstable (imaginary) phonon branches in the DFT and MEAM calculations. 
The imaginary phonon for T-[110] at the N-point corresponds to the /3 to a transformation mechanism. The imaginary L-|[lll] 
phonon is responsible for the /3 to uj transformation. 



Molecular dynamics simulations for the defects confirm 
their stability in MEAM. Calculations for larger simula- 
tion cells with 1080 and 1296 atoms show that the resid- 
ual finite-size error of the defect formation energies is 
smaller than 0.1 eV. Some interstitials can lower their 
energy in MEAM by symmetry breaking. In the a phase, 
the octahedral interstitial moves to an off-center octahe- 
dral site with an energy 0.5 eV lower than the central 
octahedral site. The dumbbell cants, lowering its energy 
by 0.09 eV. In u> the octahedral and hexahedral intersti- 
tials reduce their energy by about 1 and 0.9 eV, respec- 
tively, moving into an off-center position. No attempt 
is made to determine the defect stability in DFT or TB 
calculations. 

Experimental values for defect formation energies in 
Ti are rather difficult to obtain due to the presence of 
the a-[3 transformation and the sensitivity of diffusion 
to impurities such as oxygen. Based on an empirical re- 
lationship that connects the onset of positron trapping 
with the vacancy formation energy, positron annihilation 
experiments estimate for a-Ti a value of 1.27±0.05 eV42 
Isotope diffusion measurements result in diffusion activa- 
tion energy of 1.75 eV4^ Assuming a vacancy mechanism 
of diffusion, the two measurements lead to an estimate 
for the energy barrier of diffusion of about 0.5 eV. Both 
experimental values are significantly smaller than DFT 
predictions. The origin of the discrepancy remains un- 
clear and beyond the scope of this paper. 

The RMS force matching errors are 25% (/?), 27% 
(a) and 27% (u>). The forces depend quadratically on 
the phonon frequencies, hence, the expected error in the 
phonons is approximately half the force-matching error 
and should be of the order of 15% for all three phases. 
Errors in the low-frequency acoustic branches are signifi- 



cantly less reflecting the accuracy of the elastic constants. 



D. Accuracy of the MEAM potential 

We test the accuracy of the MEAM potential by com- 
paring the phonon spectra, surface and stacking fault en- 
ergies and energy barriers for homogeneous martensitic 
transformations to DFT, TB and experiment. 



1. Phonons 

Figure [3] compares the phonon spectra obtained with 
the MEAM potential for the three Ti phases a, (3 and 
uj with the available experimental data and results from 
DFT calculations. For all three phases, the MEAM po- 
tential reproduces the GGA phonons within about 15%, 
good agreement that can be attributed to the force- 
matching method and fitting to elastic constants. The 
acoustic branches are better reproduced than the the op- 
tical modes, reflecting the accuracy of the elastic con- 
stants. Both acoustic and optical phonons are needed to 
describe the shuffle and strain degree of freedom in the 
martensitic phase transformations. 

For the a phase the MEAM potential reproduces the 
overall trend of the experimental phonon branches. The 
optical L-[001] phonon at T is 20% too high in MEAM 
and shows the wrong curvature away from T. For the K 
and M point, the MEAM potential underestimates the 
experimental phonon frequencies by about 20 to 30%. 

No experimental data is available for the high-pressure 
uj phase. We compare the phonon spectrum for the 
MEAM potential with results from DFT calculations and 
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TABLE IV: Comparison of fitted point defect energies in the 
a and lo phases for the MEAM potential with DFT and three 
tight-binding potentials.— ' 40 ' 41 The defect energies are mea- 
sured relative to the corresponding a or u phases at about 
1 at.% defect concentration. The RMS errors of defect ener- 
gies relative to GGA are given for both phases. The tetra- 
hedral interstitial and the dumbbell- [0001] oriented along the 
c-axis are close in structure such that the tetrahedral defect 
in MEAM and GGA relaxes into the dumbbell structure. For 
the tight-binding potential by Rudin et al~ these two defects 
and the tetrahedral interstitial in u) collapse, an explanation 
and fix for this problem can be found in Ref. |40| . In the NRL 
potential the tetrahedral defect is unstable and relaxes to the 
octahedral defect. 



Site 


MEAM 


GGA 


Trmkle 4 ^ 


Rudini NRl4i 




— Q 


-defects 








Vacancy 


2.24 


2.03 


1.81 


1.92 


1.51 


Divacancy-AB" 


4.00 


3.92 


3.82 


3.68 


3.73 


Octahedral 


2.64 


2.58 


2.89 


2.55 


1.31 


Tetrahedral 


dumb. 


dumb. 


2.86 


coll. 


octa. 


Dumbbell- [0001]" 


2.21 


2.87 


2.81 


coll. 


1.81 


RMS deviation 


0.24 




0.14 


0.16 


0.76 




— LO 


-defects 








Vacancy A 


2.78 


2.92 


2.85 


3.25 


2.99 


Vacancy B 


0.82 


1.57 


1.34 


1.90 


1.01 


Octahedral 


3.79 


3.76 


4.11 


3.67 


3.20 


Tetrahedral 


2.93 


3.50 


3.58 


coll. 


2.86 


Hexahedral" 


3.31 


3.49 


3.86 


4.37 


3.20 


RMS deviation 


0.42 




0.25 


0.5 


0.42 



"Not included in potential fitting. 



find agreement between the MEAM potential and the 
DFT results comparable to that of the a phase with 
closely matching acoustic modes and larger deviations 
for the optical branches. In contrast to the a phase, for 
the lo phase the deviations for the optical branches are 
smaller and more uniform across the Brilloin zone. 

The high-temperature P phase becomes mechanically 
unstable at lower temperatures and shows a soft mode 
in the experimental data for the L-|[lll] phonon. The 
zero-temperature phonon results for the MEAM poten- 
tial and DFT reflect this instability in an unstable (imag- 
inary) phonon branch. This mode is responsible for the 
(111) plane collapse mechanism of the f3 to lo transfor- 
mation. In addition, MEAM and DFT show an unstable 
phonon branch in the T-[110] direction at the N-point 
which corresponds to the Burgers mechanism of the (3 to 
a transformation. 



2. Surface and stacking fault energies 

Large changes of coordination number provide a chal- 
lenging test for atomistic potentials. Specifically relevant 
to experiment are tests on free surfaces. Relaxations of 
low-index surfaces of the a and lo phases with MEAM 



TABLE V: Surface energies of the a and lo phase in MEAM 
compared to DFT calculations. 







a 










a (meV/A 2 ) 


[1120] 


[1100] 


[0001] 


[1120] 


[1100] 


[0001] 


MEAM 


105 


97 


92 


98 


115 


116 


GGA 


117 


153 


121 


152 


136 


133 



and DFT determine the accuracy of the classical poten- 
tial here: Calculations for increasingly larger slabs show 
that a slab thickness of more than 15 A results in surface 
energies accurate to about 1 me v/A 2 . 

Table [V] compares the surface energies of the a and 
u> phases in MEAM with DFT calculations. The overall 
agreement of the MEAM surface energies with the DFT 
values for both phases is quite remarkable considering 
the fact that free surfaces were not used to optimize the 
potential. The average MEAM surface energy is about 
20% too small. The good agreement encourages the po- 
tential's application to model systems with free surfaces 
such as voids or cracks. 

Stacking faults in hep materials alter the structural 
sequence of atomic planes in the c-direction. Their en- 
ergies test the accuracy of the potential under changes 
of bond direction and second nearest neighbor coordina- 
tion. There are three basic stacking faults possible in 
hep materials. Intrinsic stacking faults I\ and I2 change 
the hep stacking sequence ABAB to ABAB\CBCB and 
ABAB\CACA, respectively. Extrinsic stacking faults in- 
troduce additional layers in the hep stacking sequence 
such as ABAB\C\ABAB. The intrinsic stacking fault 
I2 describes a crystal sheared by a partial lattice vector 
while both the intrinsic stacking fault Ii and the extrin- 
sic stacking fault require a diffusive process. In hep ma- 
terials the I2 stacking fault energy determines the dis- 
sociation of dislocation on the basal plane into partial 
dislocation s 4 ^ 46 and has been measured for a-Ti4£ 

The MEAM potential correctly predicts a metastable 
I2 stacking fault with a high stacking fault energy 
(170 mJ/m 2 ), though not as high as in DFT (320 mJ/m 2 ) 
and experiment (300 mJ/m 2 )^ The high stacking fault 
energy in MEAM would result in a narrow splitting of 
dislocations. Elasticity predicts a splitting on the order 
of 7 A (two Burgers vectors) for a basal dislocation; this 
is consistent with a prediction of prismatic slip, as ex- 
pected for Tij 48 i 49 Since the MEAM potential predicts a 
small dislocation splitting and reproduces the elastic con- 
stants, the potential may provide an accurate description 
of dislocation interactions at short and long distances. 



3. Energy barrier from a to lo 

Figure [4] compares the energy barriers for different 
mechanism of the a to lo transformation in MEAM and 
the TB potential of Trinkle et The energy barrier 
is calculated for a two-dimensional reduced phase space 
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TB landscape barrier 

FIG. 4: (color online) Comparison of energy barriers of differ- 
ent pathways for the a to lo transformatio n 15 ' 16 with MEAM 
and TB.= The inset shows the distribution of the deviations 
between the TB and the MEAM potential. The RMS devia- 
tion is 5.5% and the maximum deviation 27% over the set of 
977 pathways. 

for all the pathways considered in Refs. [TH and[l6|. The 
first degree of freedom describes the strain of the a cell 
into the u> cell and the second degree of freedom describes 
the shuffle motion of all atoms from their a to lo position 
within the cell. The energy barriers predicted by MEAM 
agree within 27% with the highly accurate TB potential 
by Trinkle et a/4°- for all pathways and show an RMS de- 
viation of only 5.5%. This excellent agreement indicates 
a highly accurate representation of the energy landscapes 
by the MEAM potential. 



III. APPLICATION TO PHASES AND PHASE 
TRANSFORMATIONS 

The demonstrated high accuracy of the MEAM po- 
tential and its computational efficiency enable medium 
to large scale predictive simulations for titanium. This 
section presents molecular dynamics simulations for the 
different Ti phases and for an interface between the a 
and lo phases. The simulations determine the phase sta- 
bility and pressure-temperature phase diagram of a, (3 
and lo Ti and the interfacial energy and mobility for the 
a-Lo martensitic transformation. 



A. Equilibrium phase diagram of titanium 

Molecular dynamics simulations determine the stabil- 
ity of the Ti phases as a function of pressure and temper- 
ature. The simulations use the TPN (constant temper- 




Pressure [GPa] 



FIG. 5: (color online) Equilibrium phase diagram of Ti as a 
function of pressure and temperature. The MEAM potential 
accurately describes all three solid phases, a, /3 and lo, of Ti 
and captures the martensitic phase transformations between 
them. Near the triple point, the a, f3 and to phases are nearly 
degenerate. As a result, whether the /3 phase transforms to 
either a or to is controlled by the initial conditions of the 
molecular dynamics simulations. 

ature, constant pressure, constant number) ensemble 50 
with a time step of 1 fs. To estimate the stability range 
of the a, [3 and to phases we perform molecular dynamics 
starting from a cubic cell of (3 with 432 atoms, that is 
comensurate with all three phases if properly strained. 
For each pressure and temperature value we simulate 
up to 1 ns and observe the phase evolution of the sys- 
tem. Simulations for a solid-liquid interface containing 
864 atoms estimate the melting temperature of the j3 
phase. 

Figure O shows the predicted Ti equilibrium phase di- 
agram as a function of pressure and temperature for the 
classical MEAM potential. At pressures below 7 GPa 
and temperatures below 1200 K the j3 phase transforms 
into the a phase by a shear and shuffle motion of the 
atoms. At pressures above 8 GPa and temperatures be- 
low about 1300 K the (3 phase transforms into the lo 
phase. The transition temperature between the a and 
f3 phase is nearly independent of the pressure while the 
transition temperature for the (3 to lo transition increases 
with pressure. The triple point between the a, (3 and lo 
phase occurs at about 8 GPa and 1200 K. At zero pres- 
sure the f3 phase melts near 1900 K. The melting tem- 
perature first increases with pressure up to about 2000 K 
at 4 GPa and then slowly decreases with pressure. 

The phase diagram of the MEAM potential agrees 
closely with experimental observations. The a- (3 transi- 
tion occurs in experiment at 1155 K compared to 1250 K 
in the simulation. The (3 phase melts at 1943K in close 
agreement with the MEAM value of 1900 K. Measure- 
ments of the a-u transformation pressure show a large 
hysteresis with a transformation onset ranging from 2.9 
to 9.0 GPa£ i 51 ' 52 i 53 The accepted equilibrium transfor- 
mation pressure of 2.0 ± 0.3 GPa was estimated from 
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FIG. 6: (color online) Molecular dynamics simulations of the a-uj martensitic transformations, (a) The relaxed a-uj interface 
is created using the TAO-1 pathway and relaxed to produce the initial interfaces. The energy for the interface is 30 meV/A 2 . 
Unidentified atoms are colored gray, while a and u) atoms are colored blue and green, respectively, (b) Transformed geometry 
at zero compression. Molecular dynamics simulations of the relaxed interface structure at 1200 K for 200 ps transform the 
entire to phase region. When the structure is fully relaxed to zero temperature, it is identified as entirely a, except for the 
small region where the two interfaces merge and leave behind a number of defects, (c) Transformed geometry at 10% volume 
compression. The relaxed interface structure is first compressed by 10% (corresponding to ~10 GPa) and then simulated at 
1200 K for 200 ps. The entire a region transforms to lu, except for a number of point defects (interstitials and vacancies). For 
illustration, the slab is rotated around the horizontal axis so that the c axis of lu is perpendicular to the page. 



samples under shear stress that reduce the hysteresis^ 
Experimental values for the triple point range from 8 GPa 
to 9 GPa and 900 K to 1100 K^i^ similar to the MEAM 
potential values of 8 GPa and 1200 K. The close agree- 
ment of the MEAM phase diagram with experimental 
data enables quantitative simulations for the martensitic 
phase transformations between the a, (3 and lu phases. 



B. Pressure-induced martensitic phase 
transformations 

Under pressure Ti transforms from a to lu via the TAO- 
1 mechanismi 15 i 16 As a first step towards a detailed un- 
derstanding of the nucleation, we simulate a-Lu interfaces 
under compression at finite temperature. The growth of 
a nucleus of a daughter phase in a parent phase is con- 
trolled by the mobility of the interface at finite temper- 
ature under a driving force. For the Ti a to lu marten- 
sitic phase transformation the free enthalpy difference be- 
tween the two phases provides the driving force. 

We construct a-ou interfaces using the TAO-1 supercell 
to study the dynamics of the martensitic phase trans- 
formation. We set up interfaces between periodic slabs 
of untransformed TAO-1 supercells (a phase) and trans- 
formed TAO-1 supercells {lu phase) that minimize lattice 
mismatch and strain while retaining periodicity. The re- 



sulting interfaces are consistent with the pathway, mini- 
mize mismatch at the boundaries and minimize strain in 
each phase. The simulation cell has periodic boundary 
conditions in all three dimensions and contains a total of 
3,600 atoms, half in each of the two phases. The system 
consists of alternating a and lu layers of 100 A thickness. 
Relaxations of the initial interface estimate an interfacial 
energy between a and lu of 30 mc v/A 2 , roughly a third 
of the calculated surface energies (see Tab. N} . 

All simulations are performed with OttMMS^i using a 
frozen cell geometry and a Langevin thermostat to pro- 
duce a constant temperature of 1200 K. A time step of 
1 fs is used for all numerical integration with the velocity- 
Verlet propagator. 

Figure [6] shows the results of the molecular dynam- 
ics simulations for this cell at 1200 K. The simulations 
are performed at two different volumes: zero compres- 
sion for the Lu—ta transformation corresponding to ap- 
proximately GPa and 15% volume compression for the 
a^LU transformation corresponding to about 15 GPa. At 
both pressures the interface between a and lu is mobile. 
At GPa the system transforms completely to a and at 
15 GPa completely to lu, both within only 200 ps. In 
both cases the interfaces between a and lu approach each 
other and partially annihilate, leaving behind a number 
of interstitial and vacancy defects. 

Temperature alone does not drive the transformation. 
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We performed runs using the a and the u> slabs by them- 
selves at 1200 K with no compression and 10% compres- 
sion for 1 ns. In both cases the initial structure remained 
for the duration of the simulation. This indicates that ho- 
mogeneous nuclcation is not likely on the time scale of 
nanoseconds in cells with a few thousand atoms, while 
the motion of the interface does occur on such short time 
scale in these cells. 



IV. CONCLUSION 

We developed and tested a classical potential for the 
complex phase transformations of the technologically im- 
portant Ti system. The potential is of the modified 
embedded atom form ensuring computational efficiency 
with parameters optimized to density functional calcu- 
lations. The optimized potential describes the structure 
and energetics of all three phases of Ti, the a, (3 and u> 
phases. The elastic constants, phonon frequencies, sur- 



face energies and defect formation energies closely match 
density functional results even when these were not in- 
cluded in the fitting procedure. 

Molecular dynamics simulations of the phase stability 
determine the potential's equilibrium phase diagram in 
close agreement with experimental measurements. Simu- 
lations for the mobility of an a-u> interface demonstrate a 
high intcrfacial mobility corresponding to the martensitic 
character of the a—u> transformation. The potential en- 
ables quantitative studies of point defect evolution, grain 
boundary structures and mobility, as well as phase trans- 
formations in the Ti system. 
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